% Solves the model in the paper. Takes some parameter from outside, takes some
% additional moments from the data and chooses the remaining parameters to
% match these moments. Then conducts counterfactual exercises across Indian
% states, India over time, Indis Vs US, and 3 richest and 3 poorest states

clear all; close all;

%% Parameters and Targets

min_q_vec = [1;0.1];
lov_vec = [0;1]
% min_q_vec = [1];
% lov_vec = [0]

for min_q_loop = 1:size( min_q_vec,1)
for lov_loop = 1:size(lov_vec,1)

close all;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters set externally/chosen from literature %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

min_q = min_q_vec(min_q_loop); % The 'q' of the worst quality. ALLOWING FOR DIFFERENT VALUES
sigma = 5;              % Elasticity of substitution between different varities of the same quality
h = .24;                % Share of population that is skilled
N = 400000000;          % Population - Normalization
sigma_us = 1.75*1;      % Elasticity of Substitution between high and low skilled workers
delta = 0.1;            % Exogenous exit rate for firms - ONLY THERE FOR GENERALITY. NOT BEING USED AS INTERPRETING THINGS AS STATIC
beta = 1;               % Rate of time discounting - ONLY THERE FOR GENERALITY. NOT BEING USED AS INTERPRETING THINGS AS STATIC
mean_entrant = 0;       % Mean of normal from which log(A) is drawn. Should be a normalization.
A_S = 1;                % Skilled Productivity - ONLY THERE FOR GENERALITY. NOT BEING USED
z_input_q = 0;          % Amount of homogenous goods needed for each quality - ONLY THERE FOR GENERALITY. NOT BEING USED
A_Z = 1;                % Neutral productivity of homogenous good - ONLY THERE FOR GENERALITY. NOT BEING USED
theta_Z = 0.5;          % Share of unskilled in homogenous production - ONLY THERE FOR GENERALITY. NOT BEING USED
% Love of Variety Parameter - Robustness to no Love of Variety and Full
% Love of Variety
if lov_vec(lov_loop) == 0
    sheet_name = strcat('q=',num2str(min_q),';No LoV');  % This isn't a parameter - just a string which gets used in the writing results to excel.
    eta = 1/(sigma-1);                                   % Extent of LoV - 0 means full LoV and 1/(sigma-1) means no LoV
elseif lov_vec(lov_loop) == 1
    sheet_name = strcat('q=',num2str(min_q),';LoV');
    eta = 0;                                             % Extent of LoV - 0 means full LoV and 1/(sigma-1) means no LoV
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Targets, Grid Size, and Other Parameters %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

prod_grid_size = 100;               % Number of grid points for productivity.
grid_dev_sd = 5;                    % The grid for producitivity will be +- these many standard deviations
file_name = 'results_postJM';
non_equi_space = 0;                 % If 0, then equispace
log_base = exp(1);                  % If Delta's are not equispaced, then which log base to take
w_U = 1;                            % Wage of unskilled - numeraire
w_S_target = 1.6*w_U;               % Target wage premium
price_size_slope_target = .1;       % Price ratio target based on size-price regressions from ASISUM.
price_income_slope_target = .1;     % Target for the price income relation taken from NSS household surveys.
std_dev_target = .64;               % Target for std deviation of ln L
w_S_max = 5;                        % Bisection for w_S done between 0 and this number. Increase if w_S hits upper limit 
tol = 1e-6;                         % Tolerance
% Number of quality types is chosen. In addition, choose the TARGET size
% for each quality type. I choose the target size for the worst and best
% quality such that each quality type's average size is 2 times the the
% previous qualities average size. 
no_q = 12;                     % How many quality level to allow for
min_size_target = 1.25;        % Choose 1.25 because as we double this, we get 2.5,5,10 etc
max_size_target = min(min_size_target*2^(no_q - 1),10000);     % Keep doubling for each quality level
size_q_target = exp(linspace(log(min_size_target),log(max_size_target),no_q))';    % TARGET size for each quality level
% Ratio of skilled to unskilled workers in different qualities based on
% sizes. This is based roughly on extrapolation from the NSS. Although, I
% am taking the levels and slope with size of this ratio from the NSS, what
% I end up matching in the data is ratio of these to the lowest quality
% (smallest)
U_S_2pt5 = 5;                  % Taken roughly from the NSS Schedule 10
U_S_12pt5 = 2.85;              % Taken roughly from the NSS Schedule 10
min_U_to_S_ratio = .5;
slope_U_to_S_ratio = (U_S_2pt5 - U_S_12pt5)/log(2.5/12.5); 
L_U_by_L_S_q_target = U_S_2pt5 + slope_U_to_S_ratio*log(size_q_target/2.5);
L_U_by_L_S_q_target = max(L_U_by_L_S_q_target,min_U_to_S_ratio);
% Data distogram points computed. Share of empployment in different size
% categories in the data which will be targets in the calibration.
% Target employment share for each quality level will be chosen to roughly
% match this data histogram.
sd_cdf_data = csvread('data_inputs\sd_cdf_asisum_0.csv');
size_hist_cutoffs = [5;10;20;50;100;250;500;1000;2500;5000];
for ii = 1:size(size_hist_cutoffs,1)

cdf_hist_pts_data(ii,1) = sd_cdf_data(find(sd_cdf_data(:,1)>size_hist_cutoffs(ii),1)-1,2);

end
pdf_hist_pts_target = [cdf_hist_pts_data(1,1);diff(cdf_hist_pts_data); 1- cdf_hist_pts_data(end,1) ]; % Includes the last category i.e. >5000
pdf_hist_pts_target = pdf_hist_pts_target./sum(pdf_hist_pts_target);
% Computing a guess for share of employment in different quality levels
% which will roughly match the actual size distribution. Will do a fmincon
% starting from this guess to match the histogram computed above. Note that
% the dimension of sh_emp_q_guess is 1 less than no_q, simply because the
% last sh_emp_q_guess must be 1 minus the sum of all the other
% sh_emp_q_guess.
size_target_midpts_q = (size_q_target(1:end-1) + size_q_target(2:end))/2;%size_target_q(1:end-1); % 
sd_cdf_data_interp = interp1(log(sd_cdf_data(:,1)),sd_cdf_data(:,2),log(size_target_midpts_q),'cubic');
sh_emp_q_guess(1,1) = sd_cdf_data_interp(1,1);
sh_emp_q_guess(2:no_q-1,1) = diff(sd_cdf_data_interp);




%%%%%%%%%%%%%%%
%% Functions %%
%%%%%%%%%%%%%%%

u_fn = @(q,P_q,w,a,pi)q.*log((w + pi)./P_q)+ a;  % Utility function
MC_L_fn = @(w_U,w_S,theta,A,sigma_us,A_S)1./A./...
        ((theta.^sigma_us)*(1/w_U)^(sigma_us-1) + ((1-theta).^sigma_us)*(A_S/w_S)^(sigma_us-1)).^(1/(sigma_us-1)); % Marginal Cost
L_U_fn = @(d,A,theta,w_U,w_S,sigma_us,A_S)d./A./((theta + ((A_S*w_U/w_S*(1-theta)./theta).^(sigma_us-1))...
    .*(1 - theta)).^(sigma_us/(sigma_us-1))); % Amount of unskilled labor demanded for a given level of output
L_S_fn = @(L_U,theta,w_U,w_S,sigma_us,A_S)L_U.*(((w_U/w_S*(1-theta)./theta).^sigma_us))*(A_S^(sigma_us-1)); % Skilled labor given unskilled
% Inverse of marginal cost for a CES between low and high skill without the
% productivity term
MC_L_part_fn = @(w_U,w_S,theta,sigma_us,A_S)...
        ((theta.^sigma_us)*(1/w_U)^(sigma_us-1) + ((1-theta).^sigma_us)*(A_S/w_S)^(sigma_us-1)).^(1/(sigma_us-1));
% Share of labor devoted to entry - analytical solution using the free
% entry conditions and total expenditure equals wages earned.
entry_sh_L_fn = @(sigma,h,w_S,w_U,entry_sh_S,pi,beta,delta)delta/(1-beta*(1-delta))*1/sigma*(h*w_S + (1-h)*w_U + pi)/(w_U*(1-entry_sh_S)+w_S*entry_sh_S);


% Things to input
inputs_struc.N=N;
inputs_struc.h=h;
inputs_struc.sigma_us=sigma_us;
inputs_struc.sigma=sigma;
inputs_struc.w_U=w_U;
inputs_struc.eta=eta;
inputs_struc.A_Z=A_Z;
inputs_struc.theta_Z=theta_Z;
inputs_struc.mean_entrant= mean_entrant;
inputs_struc.delta = delta;
inputs_struc.beta = beta;
inputs_struc.min_q = min_q;
inputs_struc.z_input_q = z_input_q;
inputs_struc.sd_cdf_data = sd_cdf_data;
inputs_struc.A_S = A_S;

inputs_struc.no_q=no_q;
inputs_struc.prod_grid_size=prod_grid_size;
inputs_struc.grid_dev_sd=grid_dev_sd;
inputs_struc.size_hist_cutoffs = size_hist_cutoffs;
inputs_struc.sh_emp_q_guess = sh_emp_q_guess;

inputs_struc.tol=tol;
inputs_struc.w_S_max=w_S_max;

inputs_struc.w_S_target=w_S_target;
inputs_struc.price_size_slope_target = price_size_slope_target;
inputs_struc.price_income_slope_target = price_income_slope_target;
inputs_struc.size_q_target = size_q_target;
inputs_struc.L_U_by_L_S_q_target = L_U_by_L_S_q_target;
inputs_struc.std_dev_target = std_dev_target;
inputs_struc.pdf_hist_pts_target = pdf_hist_pts_target;

inputs_struc.non_equi_space = non_equi_space;
inputs_struc.log_base = log_base;

fns_struc.u_fn = u_fn;
fns_struc.MC_L_fn = MC_L_fn;
fns_struc.L_U_fn = L_U_fn;
fns_struc.L_S_fn = L_S_fn;
fns_struc.MC_L_part_fn = MC_L_part_fn;
fns_struc.entry_sh_L_fn = entry_sh_L_fn;



%%%%%%%%%%%%%%%%%
%% Calibrating %%
%%%%%%%%%%%%%%%%%

% Bisection for sd_entrants
sd_iter = 0;
sd_min = 0.05;
sd_max = 0.15;
sd_error = 1;



while sd_error>1e-5 && sd_iter<30
    
    sd_iter = sd_iter+1;
    sd_entrant = (sd_min + sd_max)/2;
    
    % Takes as input the std dev of productivity for entrants. Given this
    % std_dev, finds the sh of employment for each quality levels that best
    % matched hist (size distribution) and then finds the std dev of
    % overall log(Size). 
    [sd_calib,sh_emp_q_target] = sd_sh_emp_q_calib_fn(inputs_struc,sd_entrant);
    
   
    sd_error = abs(sd_calib/std_dev_target-1);
    
    if sd_calib>std_dev_target
        sd_max = sd_entrant;
    else
        sd_min = sd_entrant;
    end
    
    
end

inputs_struc.sd_entrant = sd_entrant;
inputs_struc.sh_emp_q_target = sh_emp_q_target;



iter_theta = 0;
theta_min = 0;
theta_max = 1;
theta_error = 1;

tic
% Do a bisection over theta for lowest quality to fit the wage premia. This
% is the outer loop. All other theta's (and other parameters) are chosen to
% match other moments within the loop.
while theta_error>tol&&iter_theta<30

    iter_theta = iter_theta + 1;
    % Thets_B chosen to match the wage premium - the theta for the lowest
    % qualit basically
    theta_B = (theta_min + theta_max)/2;
    % Share parameter of unskilled in higher quality levels are chosen to
    % match ratio of unskilled to skilled for different quality levels 
    theta_q = 1./(1 + ((L_U_by_L_S_q_target(1,1)./L_U_by_L_S_q_target(2:end,1)).^(1/sigma_us)).*((1-theta_B)/theta_B));   
    theta_q = [theta_B;theta_q];
    theta_q_i = repmat(theta_q,1,prod_grid_size);       % Expand to grid 
    S_to_U_ratio_q = ((w_U/w_S_target*(1-theta_q)./theta_q).^(sigma_us)).*(A_S.^(sigma_us-1));  % Ratio of skilled to unskilled given the guessed theta's and the target for w_S
    entry_sh_S_q = 1 - 1./(S_to_U_ratio_q + 1);     % The share of skilled used in entry for different qualities is the same as in production
    
    p_z = 0;            % Price of homogenous good - ONLY THERE FOR GENERALITY. NOT BEING USED
    z_input_q_i = 0;    % Intermediate inputs needed - ONLY THERE FOR GENERALITY. NOT BEING USED
    
    
    % Adding guess/calibrated values of theta's and others to input structure
    inputs_struc.theta_q_i = theta_q_i;
    inputs_struc.entry_sh_S_q =     entry_sh_S_q;
    inputs_struc.theta_q = theta_q;
    inputs_struc.p_z = p_z;
    inputs_struc.z_input_q_i = z_input_q_i;
    

    % Given calibrated value of std dev, and thetas get productivty grids
    [A_q_i,pdf_A_q_i,A_scale_q,mean_test_perdiff,sd_test_perdiff] = A_q_i_grid_dist(inputs_struc,fns_struc);
    mean_test_perdiff
    sd_test_perdiff
    inputs_struc.A_q_i = A_q_i;
    inputs_struc.pdf_A_q_i = pdf_A_q_i;
   

    pi_guess = 0;       % With no dynamics (discounting) this will be 0. Haven't done calibration for the case with discounting??
    % Guess for sh of L that will get devoted to entry. The function
    % entry_sh_L_fn gives the exact answer for this sh when all qualities
    % use the same amount of skill for entry. It also gives the exact
    % answer when entry requirement if skill is exactly the same as
    % production requirement (and also setting entry skill requirement to
    % 'h'). The latter is the case in the calibration but for generality I
    % check whether this conjecture is true.
    entry_sh_L_guess = entry_sh_L_fn(sigma,h,w_S_target,w_U,h,pi_guess,beta,delta);
    error_entry_sh_L = 1;
    iter_entry_sh_L = 1;
    
   % This loop takes the targets for avg size and employment share for
   % different quality levels and backs out targets for mass of entrants
   % and then figures out fixed costs to hit the target for M_q. It then
   % checks to make sure that given the fixed costs calibrated, the guess
   % for share of labor devoted to entry is realized.
   while error_entry_sh_L>tol&&iter_entry_sh_L<30

       iter_entry_sh_L = iter_entry_sh_L + 1;
        % Some useful constants
        L_q_constant = ((w_U/w_S_target*(1-theta_q)./theta_q).^sigma_us)*(A_S^(sigma_us-1));
        MC_L_part = MC_L_part_fn(w_U,w_S_target,theta_q,sigma_us,A_S);
        M_q_target = sh_emp_q_target.*(1-entry_sh_L_guess).*N./size_q_target;           % The mass of firms we need in equilibrium to match the targeted size of each q
        L_U_q_target = sh_emp_q_target.*(1-entry_sh_L_guess).*N./(1 + L_q_constant);    % How much unskilled labor will be employed by each quality
        % The next two lines basically implement equations labelled "Entry cost
        % calib" in the file "FGH_CES_lnC_HetrProd_LoVParameter_calib.lyx"
        x_part_q = (theta_q + (1-theta_q).*(((1-theta_q)./theta_q*w_U/w_S_target*A_S).^(sigma_us-1))).^(sigma_us/(sigma_us-1));
        f_q = 1/(1-beta*(1-delta))*1/(sigma-1)*x_part_q.*(1./MC_L_part).*L_U_q_target./M_q_target./((1-entry_sh_S_q).*w_U + entry_sh_S_q.*w_S_target);

        entry_sh_L_new = sum(f_q.*M_q_target*delta)/N;
        error_entry_sh_L = abs(entry_sh_L_new/entry_sh_L_guess - 1);
        entry_sh_L_guess = entry_sh_L_new*0.5 + 0.5*entry_sh_L_guess;
    
   end
   
   % Adding guess/calibrates values of entry costs and emp share in entry.
    inputs_struc.entry_sh_L =     entry_sh_L_guess;
    inputs_struc.pi_target =     0;
    inputs_struc.f_q =  f_q;
    
    
    % Calibrate the entrants standard deviation of productivity and the max q
    % to match std dev of labor and non-homotheticity. Within this
    % calibration, all the other parameters get calibrated - this is
    % repeated right after the q_sd_solve command.
    max_q_guess = min_q + 1.5 ;
    q_sd_guess = [max_q_guess];
    q_sd_solve = q_sd_calib_fn(q_sd_guess,inputs_struc,fns_struc);


    % Vector of q using calibrated max_q
    max_q_solve = q_sd_solve(1);
    q = linspace(min_q,max_q_solve,no_q)';
    if non_equi_space==1
        q =log(linspace(log_base^(min_q),log_base^(max_q_solve),no_q))'./log(log_base);
    end  
    inputs_struc.q = q;

    % a_q, given q to match size distribution
    a_q = a_q_calib_fn(inputs_struc,fns_struc);
    inputs_struc.a_q = [0;a_q];

    % Solve for equilibrium ans make sure all targets are met!
    Output_struc = solve_equi(inputs_struc,fns_struc);
    w_S_calib=Output_struc.w_S ;

   theta_error = abs(w_S_calib/w_S_target - 1);
   
   if w_S_calib>w_S_target
       theta_min = theta_B;
   else
       theta_max = theta_B;
   end
   
   
    
end

toc

Output_struc = solve_equi(inputs_struc,fns_struc);

% Stuff which will be written into excel
write_excel = [inputs_struc.min_q;inputs_struc.h;1;1;Output_struc.w_S;Output_struc.L_U_clearing;Output_struc.L_S_clearing];
rel_price_write = Output_struc.P_q_norm;
rel_price_write = rel_price_write./rel_price_write(1);
avg_size_write = (Output_struc.L_U_q + Output_struc.L_S_q)./Output_struc.M_q;
theta_write = inputs_struc.theta_q_i;
theta_write = theta_write(:,1);
rel_prod_write = inputs_struc.A_q_i;
rel_prod_write = rel_prod_write(:,1)./rel_prod_write(1,1);
sh_S_entry_write = inputs_struc.entry_sh_S_q;
write_excel = [write_excel;0;0;rel_price_write; avg_size_write; theta_write; rel_prod_write ;sh_S_entry_write];
Output_struc.write_excel = write_excel;


w_S_calib=Output_struc.w_S ; L_U_clearing_calib=Output_struc.L_U_clearing; L_S_clearing_calib=Output_struc.L_S_clearing;
p_q_i_calib=Output_struc.p_q_i;MC_L_q_i_calib=Output_struc.MC_L_q_i; P_q_calib=Output_struc.P_q; 
rho_q_U_calib=Output_struc.rho_q_U;rho_q_S_calib=Output_struc.rho_q_S; D_q_calib=Output_struc.D_q; 
M_q_calib=Output_struc.M_q; P_q_norm_calib = Output_struc.P_q_norm;
d_q_i_calib=Output_struc.d_q_i; l_U_q_i_calib=Output_struc.l_U_q_i; l_S_q_i_calib=Output_struc.l_S_q_i;
L_U_q_calib=Output_struc.L_U_q; L_S_q_calib=Output_struc.L_S_q; L_entry_U_q_calib=Output_struc.L_entry_U_q;
L_entry_S_q_calib=Output_struc.L_entry_S_q; free_entry_q_calib=Output_struc.free_entry_q;
n_q_calib = Output_struc.n_q; pi_calib = Output_struc.pi;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Check that Calibration targets are met %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% The slope on the consumer side i.e. how much higher an average price are
% high skilled paying. Try different weighting schemes and allow for direct
% consumption of intermediate goods by the consumer. When allowing for
% direct consumption of intermediate goods, does not matter whether I do
% expenditure weighting or not. When consumping the final good, then the
% slope is the same when using the normalized prize (without LoV).
exp_wt = d_q_i_calib.*p_q_i_calib./repmat(sum(d_q_i_calib.*p_q_i_calib,2),1,prod_grid_size);
dd_wt_S = repmat(rho_q_S_calib,1,prod_grid_size)./sum(sum(repmat(rho_q_S_calib,1,prod_grid_size)));
dd_wt_U = repmat(rho_q_U_calib,1,prod_grid_size)./sum(sum(repmat(rho_q_U_calib,1,prod_grid_size)));
slope_no_wt  = (sum(sum(dd_wt_S.*log(p_q_i_calib))) - sum(sum(dd_wt_U.*log(p_q_i_calib))))/(log(w_S_calib) - log(w_U));
exp_wt_S = dd_wt_S.*exp_wt/sum(sum(dd_wt_S.*exp_wt));
exp_wt_U = dd_wt_U.*exp_wt/sum(sum(dd_wt_U.*exp_wt));
slope_exp_wt = (sum(sum(exp_wt_S.*log(p_q_i_calib))) - sum(sum(exp_wt_U.*log(p_q_i_calib))))/(log(w_S_calib) - log(w_U));
slope_exp_wt = (sum(sum(repmat(rho_q_S_calib,1,prod_grid_size).*log(p_q_i_calib).*exp_wt)) - sum(sum(repmat(rho_q_U_calib,1,prod_grid_size).*log(p_q_i_calib).*exp_wt)))/(log(w_S_calib) - log(w_U));
slope_final_good = (sum(rho_q_S_calib.*log(P_q_norm_calib)) - sum(rho_q_U_calib.*log(P_q_norm_calib)))/(log(w_S_calib) - log(w_U));

% Std deviation of employment
size_q_i_calib  = l_U_q_i_calib + l_S_q_i_calib;
ln_size_1dim_calib = log(reshape(size_q_i_calib(:,:),[no_q*prod_grid_size,1]));
mass_firms_q_i_calib = pdf_A_q_i.*repmat(M_q_calib,[1,prod_grid_size]);
mass_firms_1dim_calib = reshape(mass_firms_q_i_calib,[no_q*prod_grid_size,1]);
avg_ln_size_calib = sum((ln_size_1dim_calib.*mass_firms_1dim_calib))/sum(mass_firms_1dim_calib);
std_dev_ln_size_calib = sqrt(sum(((ln_size_1dim_calib.^2).*mass_firms_1dim_calib))/sum(mass_firms_1dim_calib) - avg_ln_size_calib^2);
error_sd = abs(std_dev_ln_size_calib/std_dev_target - 1)

% Wage Premium
dev_w_S_calib = w_S_calib/w_S_target - 1

% Ratio of skilled to unskilled for different size categories
L_U_by_L_S_reltosmall_q_calib = (L_U_q_calib./L_S_q_calib)./(L_U_q_calib(1,1)./L_S_q_calib(1,1));
dev_L_U_by_L_S_reltosmall_q_calib = max(abs(L_U_by_L_S_reltosmall_q_calib./(L_U_by_L_S_q_target./L_U_by_L_S_q_target(1,1))-1))

% Price ratio across qualities - CHECK
price_ratio_target_q  = exp(price_size_slope_target*log(size_q_target));  
price_ratio_target_q = price_ratio_target_q/price_ratio_target_q(1,1);
dev_avg_price_q_calib = max(abs(P_q_norm_calib./price_ratio_target_q - 1))

% Average size across qualities
dev_avg_size_q_calib = max(abs(((L_U_q_calib + L_S_q_calib)./M_q_calib)./(size_q_target) - 1))

% Sh of employment in different quality levels - eventually helps matching
% size distribution
sh_L_prod_q_calib = (L_U_q_calib + L_S_q_calib)/(sum(L_U_q_calib + L_S_q_calib));
dev_sh_L_prod_q_calib = max(abs(sh_L_prod_q_calib./sh_emp_q_target-1))

% Test that variance of employment within size category matches the
% theoretical variance given sigma and entrants productivity variance.
ln_size_q_i_calib = log(size_q_i_calib);
avg_ln_size_calib_q = sum(ln_size_q_i_calib.*mass_firms_q_i_calib,2)./sum(mass_firms_q_i_calib,2);
var_ln_size_calib_q = sum((ln_size_q_i_calib.^2).*mass_firms_q_i_calib,2)./sum(mass_firms_q_i_calib,2) - avg_ln_size_calib_q.^2;
test_var_ln_size_calib_q = var_ln_size_calib_q/(((sigma-1)^2)*(sd_entrant^2))-1;

Output_struc.pcgdp_calib_rel = 1;
Output_struc_calib = Output_struc;


%%%%%%%%%%%%%%%%
%% Statistics %%
%%%%%%%%%%%%%%%%

size_hist_cutoffs_graph = [5;10;20;50;100;250;500;1000;2500];
inputs_struc.size_hist_cutoffs_graph = size_hist_cutoffs_graph;
size_1dim_calib  = exp(ln_size_1dim_calib);
for i = 1:size(size_hist_cutoffs_graph,1)-1
    
    indicator_model = (size_1dim_calib<=size_hist_cutoffs_graph(i));
    indicator_data = find(sd_cdf_data(:,1)>=size_hist_cutoffs_graph(i),1,'first');
    size_hist_cdf_model_calib(i,1) = sum(indicator_model.*size_1dim_calib.*mass_firms_1dim_calib)/sum(size_1dim_calib.*mass_firms_1dim_calib);
    size_hist_cdf_data(i,1) = sd_cdf_data(indicator_data,2);
    
end
size_hist_cdf_model_calib(i+1,1) = 1;
size_hist_cdf_data(i+1,1) = 1;

size_hist_model_calib = [size_hist_cdf_model_calib(1,1);diff(size_hist_cdf_model_calib)];
size_hist_data = [size_hist_cdf_data(1,1);diff(size_hist_cdf_data)];
 
current_folder = cd;

if min_q_loop==1&&lov_loop==1
sd_model_fit = figure(1);
sd_model_fit_sr1 = bar([size_hist_data,size_hist_model_calib]);
sd_model_fit_ch = get(sd_model_fit_sr1,'Children');
set(sd_model_fit_ch{1},'FaceColor',[186/255 200/255 212/255]);
set(sd_model_fit_ch{2},'FaceColor',[144/255 53/255 59/255]);
% hatchfill(sd_model_fit_sr1);
% title('Size Distribution - Data vs Model');
handle = get(gca, 'title');
set(handle, 'FontWeight', 'Bold','FontSize',12);
xlabel('Size');
handle = get(gca, 'xlabel');
set(handle, 'FontWeight', 'Bold','FontSize',15,'FontName','TimesNewRoman');
ylabel('Share of Employment');
handle = get(gca, 'ylabel');
set(handle, 'FontWeight', 'Bold','FontSize',15,'FontName','TimesNewRoman');
% xticklabel_rotate([1:size(size_hist_cutoffs_graph,1)],45,{'<=5', '5 to 9 ', '10 to 19',  '20 to 49', '50 to 99',...
%     '100 to 249', '250 to 499', '500 to 999', '>1000'},'interpreter','none');
% xtickangle([1:size(size_hist_cutoffs_graph,1)],45,{'<=5', '5 to 9 ', '10 to 19',  '20 to 49', '50 to 99',...
%     '100 to 249', '250 to 499', '500 to 999', '>1000'},'interpreter','none');
[legend_handle] = legend('Data','Model');
set(legend_handle,'FontName','TimesNewRoman','FontSize',15);
% legend_handle_ch = get(legend_handle,'Children');
% legend_handle_grand_ch = get(legend_handle_ch(3),'Children');
% set(legend_handle_grand_ch,'FaceColor',[186/255 200/255 212/255]);
% legend_handle_grand_ch = get(legend_handle_ch(1),'Children');
% set(legend_handle_grand_ch,'FaceColor',[144/255 53/255 59/255]);
print(sd_model_fit,'-depsc2',strcat(current_folder,'\results\model_fit.eps'));
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Counterfactual - Indian States %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Bring in the sh of employment in small plants for different states. The
% columns represent: state_code, fitted sh emp<5, pcSGDP relative to India,
% india's pc GDP, and actual sh emp<5...
sd_less5_allSt = csvread('data_inputs\matlab_st_sd_relGDP.csv');
diff_data = max(sd_less5_allSt(:,2)) - min(sd_less5_allSt(:,2));        % Difference in sh of emp in <5 for poorest vs richest state based on regresseion fitted values


pcGDP_st_rel_Ind_vec = [min(sd_less5_allSt(:,3));max(sd_less5_allSt(:,3))];   % pcSGDP of richest and poorest state relative to India as a whole
pcGDP_st_vec = sd_less5_allSt(1,4)*pcGDP_st_rel_Ind_vec;                      % pcSGDP of richest and poorest state
h_hat_vec = [.13;0.43]; 	% Supply of skill of richest and poorest state
pcGDP_vec_target = pcGDP_st_rel_Ind_vec;   % pcGDP of richest and poorest state

% Computing pcGDP of calibrated baseline and storing calib prices and pcGDP
% into input structure
pcgdp_calib_own_prices = sum(sum(p_q_i_calib.*d_q_i_calib.*pdf_A_q_i.*repmat(M_q_calib,[1,prod_grid_size])))/N;    % pcGDP of calib baseline
pcgdp_calib_rel = pcgdp_calib_own_prices./pcgdp_calib_own_prices;   % pcGDP of calib baseline relative to calib baseline
inputs_struc.p_q_i_calib = p_q_i_calib;
inputs_struc.pcgdp_calib_own_prices = pcgdp_calib_own_prices;



% SOlving for the counterfactual, adjusting thetas, h and productivity to
% match supply og skill, wage premia, and pcGDP.
for i = 1:1:size(pcGDP_vec_target,1)

    % New supply of skill and target GDP in counterfactual
    inputs_struc.pcgdp_target = pcGDP_vec_target(i);
    inputs_struc.h = h_hat_vec(i);
    
    [Output_struc,inputs_struc] = solve_counterfactual(inputs_struc,fns_struc);
    
    % Make sure that errors are around 0
    d_q_i_C=Output_struc.d_q_i;
    M_q_C=Output_struc.M_q;
    pcgdp_prodvec_C_calib_prices = sum(sum(p_q_i_calib.*d_q_i_C.*pdf_A_q_i.*repmat(M_q_C,[1,size(p_q_i_calib,2)])))/N;
    pcgdp_prodvec_C_rel_calib_prices = pcgdp_prodvec_C_calib_prices/pcgdp_calib_own_prices;
    error_prod_mult = abs(pcgdp_prodvec_C_rel_calib_prices/pcGDP_vec_target(i) - 1)
    error_theta = abs(Output_struc.w_S./w_S_target - 1)
    
    Output_struc_across_st_C{i,1} = Output_struc;
    inputs_struc_across_st_C{i,1} = inputs_struc;

end


% Write C across states to excel
excel_column = 1;
prod_mult = 1;
% Write the labels on the right most column
xlswrite(file_name,{'poorest','calib','richest'},sheet_name,strcat(char('A'+excel_column),'1'));
write_label = {'q_1';'h';'prod_mult';'pcGDP';'w_S';'L_U clearing';'L_S_clearing';'ch_sd';'sh_emp_small'};
write_label{end+1,1} = 'rel_price';
write_label{end+no_q,1} = 'avg_size';
write_label{end+no_q,1} = 'theta';
write_label{end+no_q,1} = 'rel_prod';
write_label{end+no_q,1} = 'sh S entry';
xlswrite(file_name,write_label,sheet_name,'A2');
% Write output
xlswrite(file_name,[Output_struc_across_st_C{1,1}.write_excel,Output_struc_calib.write_excel,Output_struc_across_st_C{2,1}.write_excel],sheet_name,'B2');
diff_model_across_st = min(Output_struc_across_st_C{1,1}.size_hist_cdf_model_C) - min(Output_struc_across_st_C{end,1}.size_hist_cdf_model_C);
sh_exp_less5_across_st = diff_model_across_st/diff_data
xlswrite(file_name,repmat(sh_exp_less5_across_st,1,size(pcGDP_vec_target,1)+1),sheet_name,'B9');
xlswrite(file_name,[min(Output_struc_across_st_C{1,1}.size_hist_cdf_model_C),min(size_hist_cdf_model_calib),min(Output_struc_across_st_C{2,1}.size_hist_cdf_model_C)],sheet_name,'B10');

% FIGURE ACROSS STATES%
if min_q_loop==1&&lov_loop==1
plot_pcGDP = log(sd_less5_allSt(1,4)) + log([Output_struc_across_st_C{1,1}.pcgdp_calib_rel;Output_struc_calib.pcgdp_calib_rel;Output_struc_across_st_C{end,1}.pcgdp_calib_rel]);
plot_sd = [min(Output_struc_across_st_C{1,1}.size_hist_cdf_model_C);min(size_hist_cdf_model_calib); min(Output_struc_across_st_C{end,1}.size_hist_cdf_model_C)];
sd_scatter_AllSt = figure(3);
plot_data_line = plot(log(sd_less5_allSt(:,3)) + log(sd_less5_allSt(1,4)),sd_less5_allSt(:,2),'LineWidth',2,'Color','blue');
hold on
%plot_data_pts = scatter(log(sd_less5_allSt(:,3)) + log(sd_less5_allSt(1,4)),sd_less5_allSt(:,5),75,'.','MarkerEdgeColor','b');
plot_model_line = plot(plot_pcGDP ,plot_sd,'r--','LineWidth',2);
%scatter(log(sd_less5_allSt(1,4)),size_hist_model_calib(1,1),75,'r','LineWidth',2.5);
% title('India Across States - Data vs Model');
handle = get(gca, 'title');
set(handle, 'FontWeight', 'Bold','FontSize',12);
xlabel('ln(pcGDP)');
handle = get(gca, 'xlabel');
set(handle, 'FontWeight', 'Bold','FontSize',15,'FontName','TimesNewRoman');
ylabel('Share of Employment in <=5');
handle = get(gca, 'ylabel');
set(handle, 'FontWeight', 'Bold','FontSize',15,'FontName','TimesNewRoman');
handle = legend([plot_data_line plot_model_line ],{'Data','Model'});
set(handle,'FontSize',15,'FontName','TimesNewRoman');
% set(gca, 'YLim',[0.4 0.95]);
set(gca,'YTick',[0.4:.1:0.9]);
print(sd_scatter_AllSt,'-depsc2',strcat(current_folder,'\results\sd_scatter_AllSt.eps'));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Counterfactual - Richest Poorest %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Bring in the sh of employment in small plants for different states. The
% columns represent: state_code, fitted sh emp<5, pcSGDP relative to India,
% india's pc GDP, and actual sh emp<5...
sd_cdf_data_3rich = csvread('data_inputs\matlab_sd_relGDP_3rich.csv');
sd_cdf_data_3poor = csvread('data_inputs\matlab_sd_relGDP_3poor.csv');


pcGDP_st_rel_Ind_vec = [sd_cdf_data_3poor(1,3);sd_cdf_data_3rich(1,3)];;   % pcSGDP of richest and poorest state relative to India as a whole
pcGDP_st_vec = sd_less5_allSt(1,4)*pcGDP_st_rel_Ind_vec;                      % pcSGDP of richest and poorest state
h_hat_vec = [.16;0.39]; 	% Supply of skill of richest and poorest state
pcGDP_vec_target = pcGDP_st_rel_Ind_vec;   % pcGDP of richest and poorest state



% SOlving for the counterfactual, adjusting thetas, h and productivity to
% match supply og skill, wage premia, and pcGDP.
for i = 1:1:size(pcGDP_vec_target,1)

    % New supply of skill and target GDP in counterfactual
    inputs_struc.pcgdp_target = pcGDP_vec_target(i);
    inputs_struc.h = h_hat_vec(i);
    
    [Output_struc,inputs_struc] = solve_counterfactual(inputs_struc,fns_struc);
    
    % Make sure that errors are around 0
    d_q_i_C=Output_struc.d_q_i;
    M_q_C=Output_struc.M_q;
    pcgdp_prodvec_C_calib_prices = sum(sum(p_q_i_calib.*d_q_i_C.*pdf_A_q_i.*repmat(M_q_C,[1,size(p_q_i_calib,2)])))/N;
    pcgdp_prodvec_C_rel_calib_prices = pcgdp_prodvec_C_calib_prices/pcgdp_calib_own_prices;
    error_prod_mult = abs(pcgdp_prodvec_C_rel_calib_prices/pcGDP_vec_target(i) - 1)
    error_theta = abs(Output_struc.w_S./w_S_target - 1)
    
    Output_struc_RP_C{i,1} = Output_struc;
    inputs_struc_RP_C{i,1} = inputs_struc;

end


for i = 1:size(size_hist_cutoffs_graph,1)-1
    
  
    indicator_data = find(sd_cdf_data_3poor(:,1)>=size_hist_cutoffs_graph(i),1,'first');
    if isempty(indicator_data)==1
            size_hist_cdf_data_3poor(i,1) = size_hist_cdf_data_3poor(i-1,1);
    else
        size_hist_cdf_data_3poor(i,1) = sum(sd_cdf_data_3poor(1:indicator_data,2));
    end
    
    indicator_data = find(sd_cdf_data_3rich(:,1)>=size_hist_cutoffs_graph(i),1,'first');
    size_hist_cdf_data_3rich(i,1) =sum(sd_cdf_data_3rich(1:indicator_data,2));
    
end
size_hist_cdf_data_3poor(i+1,1) = 1;
size_hist_cdf_data_3rich(i+1,1) = 1;
size_hist_data_3poor = [size_hist_cdf_data_3poor(1,1);diff(size_hist_cdf_data_3poor)];
size_hist_data_3rich = [size_hist_cdf_data_3rich(1,1);diff(size_hist_cdf_data_3rich)];


% Write C across states to excel
excel_column = 1;
prod_mult = 1;
% Write the labels on the right most column
xlswrite(file_name,{'poorest','calib','richest'},strcat(sheet_name,';RP'),strcat(char('A'+excel_column),'1'));
write_label = {'q_1';'h';'prod_mult';'pcGDP';'w_S';'L_U clearing';'L_S_clearing';'ch_sd'};
write_label{end+1,1} = 'rel_price';
write_label{end+no_q,1} = 'avg_size';
write_label{end+no_q,1} = 'theta';
write_label{end+no_q,1} = 'rel_prod';
write_label{end+no_q,1} = 'sh S entry';
xlswrite(file_name,write_label,strcat(sheet_name,';RP'),'A2');
% Write output
xlswrite(file_name,[Output_struc_RP_C{1,1}.write_excel,Output_struc_calib.write_excel,Output_struc_RP_C{2,1}.write_excel],strcat(sheet_name,';RP'),'B2');
diff_model_RP = min(Output_struc_RP_C{1,1}.size_hist_cdf_model_C) - min(Output_struc_RP_C{end,1}.size_hist_cdf_model_C);
sh_exp_less5_RP = diff_model_RP/(size_hist_data_3poor(1) - size_hist_data_3rich(1))
xlswrite(file_name,repmat(sh_exp_less5_RP,1,size(pcGDP_vec_target,1)+1),strcat(sheet_name,';RP'),'B9');

size_hist_cdf_model_3poor_prodvec = Output_struc_RP_C{1,1}.size_hist_cdf_model_C';
size_hist_cdf_model_3rich_prodvec = Output_struc_RP_C{end,1}.size_hist_cdf_model_C';
size_hist_model_3poor_prodvec = [size_hist_cdf_model_3poor_prodvec(1,1);diff(size_hist_cdf_model_3poor_prodvec)];
size_hist_model_3rich_prodvec = [size_hist_cdf_model_3rich_prodvec(1,1);diff(size_hist_cdf_model_3rich_prodvec)];
size_hist_data_3poor = [size_hist_cdf_data_3poor(1,1);diff(size_hist_cdf_data_3poor)];
size_hist_data_3rich = [size_hist_cdf_data_3rich(1,1);diff(size_hist_cdf_data_3rich)];



sd_data_diff_rich_poor = size_hist_data_3poor - size_hist_data_3rich;
sd_model_diff_rich_poor = size_hist_model_3poor_prodvec - size_hist_model_3rich_prodvec;

% FIGURE Rich Minus Poor%
current_folder = cd;
if min_q_loop==1&&lov_loop==1
sd_diff_richpoor = figure(2);
sd_diff_richpoor_bar = bar([sd_data_diff_rich_poor,sd_model_diff_rich_poor]);
sd_diff_richpoor_bar_ch = get(sd_diff_richpoor_bar,'Children');
set(sd_diff_richpoor_bar_ch{1},'FaceColor',[186/255 200/255 212/255]);
set(sd_diff_richpoor_bar_ch{2},'FaceColor',[144/255 53/255 59/255]);
%title('Difference in Share of Employment: Poor - Rich');
%handle = get(gca, 'title');
xlabel('Size');
handle = get(gca, 'xlabel');
set(handle, 'FontWeight', 'Bold','FontSize',15,'FontName','TimesNewRoman');
ylabel('Share of Employment');
handle = get(gca, 'ylabel');
set(handle, 'FontWeight', 'Bold','FontSize',15,'FontName','TimesNewRoman');
xticklabel_rotate([1:size(size_hist_cutoffs_graph,1)],45,{'<=5', '5 to 9 ', '10 to 19',  '20 to 49', '50 to 99',...
    '100 to 249', '250 to 499', '500 to 999', '>1000'},'interpreter','none');
sd_diff_richpoor_legend = legend('Data','Model');
set(sd_diff_richpoor_legend,'FontSize',15,'FontName','TimesNewRoman');
sd_diff_richpoor_legend_ch = get(sd_diff_richpoor_legend,'Children');
% sd_diff_richpoor_legend_grand_ch = get(sd_diff_richpoor_legend_ch(1),'Children');
% set(sd_diff_richpoor_legend_grand_ch,'FaceColor',[144/255 53/255 59/255]);
% sd_diff_richpoor_legend_grand_ch = get(sd_diff_richpoor_legend_ch(3),'Children');
% set(sd_diff_richpoor_legend_grand_ch,'FaceColor',[186/255 200/255 212/255]);
print(sd_diff_richpoor,'-depsc2',strcat(current_folder,'\results\sd_diff_rich_poor.eps'));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Counterfactual - India Over Time %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sh <5 in dta for different years - just copy pasted in for now
sh_less5_overTime_data = [0.769;0.626;0.673;0.578];

pcGDP_st_rel_Ind_vec = [0.536113781;0.608687296;0.782702287;1.300903325];
h_hat_vec = [.14;.18;0.22;0.31]; 	% Supply of skill of richest and poorest state
pcGDP_vec_target = pcGDP_st_rel_Ind_vec;   % pcGDP of richest and poorest state


% SOlving for the counterfactual, adjusting thetas, h and productivity to
% match supply og skill, wage premia, and pcGDP.
for i = 1:1:size(pcGDP_vec_target,1)

    % New supply of skill and target GDP in counterfactual
    inputs_struc.pcgdp_target = pcGDP_vec_target(i);
    inputs_struc.h = h_hat_vec(i);
    
    [Output_struc,inputs_struc] = solve_counterfactual(inputs_struc,fns_struc);
    
    % Make sure that errors are around 0
    d_q_i_C=Output_struc.d_q_i;
    M_q_C=Output_struc.M_q;
    pcgdp_prodvec_C_calib_prices = sum(sum(p_q_i_calib.*d_q_i_C.*pdf_A_q_i.*repmat(M_q_C,[1,size(p_q_i_calib,2)])))/N;
    pcgdp_prodvec_C_rel_calib_prices = pcgdp_prodvec_C_calib_prices/pcgdp_calib_own_prices;
    error_prod_mult = abs(pcgdp_prodvec_C_rel_calib_prices/pcGDP_vec_target(i) - 1)
    error_theta = abs(Output_struc.w_S./w_S_target - 1)
    
    Output_struc_TS_C{i,1} = Output_struc;
    inputs_struc_TS_C{i,1} = inputs_struc;

end



% Write C across states to excel
excel_column = 1;
prod_mult = 1;
% Write the labels on the right most column
xlswrite(file_name,{'poorest','calib','richest'},strcat(sheet_name,';TimeSeries'),strcat(char('A'+excel_column),'1'));
write_label = {'q_1';'h';'prod_mult';'pcGDP';'w_S';'L_U clearing';'L_S_clearing';'ch_sd'};
write_label{end+1,1} = 'rel_price';
write_label{end+no_q,1} = 'avg_size';
write_label{end+no_q,1} = 'theta';
write_label{end+no_q,1} = 'rel_prod';
write_label{end+no_q,1} = 'sh S entry';
xlswrite(file_name,write_label,strcat(sheet_name,';TimeSeries'),'A2');
% Write output
xlswrite(file_name,[Output_struc_TS_C{1,1}.write_excel,Output_struc_TS_C{2,1}.write_excel,Output_struc_TS_C{3,1}.write_excel,Output_struc_calib.write_excel,Output_struc_TS_C{4,1}.write_excel],strcat(sheet_name,';TimeSeries'),'B2');
diff_model_TS = min(Output_struc_TS_C{1,1}.size_hist_cdf_model_C) - min(Output_struc_TS_C{end,1}.size_hist_cdf_model_C);
sh_exp_less5_TS = diff_model_TS/(max(sh_less5_overTime_data) - min(sh_less5_overTime_data))
xlswrite(file_name,repmat(sh_exp_less5_TS,1,size(pcGDP_vec_target,1)+1),strcat(sheet_name,';TimeSeries'),'B9');


% FIGURE Time Series %
current_folder = cd;
if min_q_loop==1&&lov_loop==1
plot_sd = [min(Output_struc_TS_C{1,1}.size_hist_cdf_model_C);min(Output_struc_TS_C{2,1}.size_hist_cdf_model_C);min(Output_struc_TS_C{3,1}.size_hist_cdf_model_C);min(size_hist_cdf_model_calib); min(Output_struc_TS_C{end,1}.size_hist_cdf_model_C)];
sd_overTime = figure(4);
bar([sh_less5_overTime_data(1:3,1);min(size_hist_cdf_data);sh_less5_overTime_data(end,1)],0.5,'FaceColor',[0.6500         0         0])
hold on
plot(plot_sd,'-o','LineWidth',3,'MarkerSize',10,'Color',[0         0    0.7500])
set(gca,'YLim',[0.4 0.8]);
% title('India Over Time - Data vs Model');
handle = get(gca, 'title');
xlabel('Year')
handle = get(gca, 'xlabel');
set(handle, 'FontWeight', 'Bold','FontSize',15,'FontName','TimesNewRoman');
ylabel('Share of Employment in  <=5')
handle = get(gca, 'ylabel');
set(handle, 'FontWeight', 'Bold','FontSize',15,'FontName','TimesNewRoman');
set(gca,'XTickLabels',{'1989','1994','2000','2005','2009'})
handle = legend('Data','Model');
set(handle,'FontSize',15,'FontName','TimesNewRoman');
print(sd_overTime,'-depsc2',strcat(current_folder,'\results\sd_india_over_time.eps'))
end
end
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Counterfactual - India vs US     %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% if min_q_loop==1&&lov_loop==1
% 
%     % Sh <5 in dta for different years - just copy pasted in for now
%     sh_less5_US_data = [0.02];
% 
%     pcGDP_st_rel_Ind_vec = [17];
%     h_hat_vec = [0.75]; 	% Supply of skill of richest and poorest state
%     pcGDP_vec_target = pcGDP_st_rel_Ind_vec;   % pcGDP of richest and poorest state
% 
% 
%     % SOlving for the counterfactual, adjusting thetas, h and productivity to
%     % match supply og skill, wage premia, and pcGDP.
%     for i = 1:1:size(pcGDP_vec_target,1)
% 
%     % New supply of skill and target GDP in counterfactual
%     inputs_struc.pcgdp_target = pcGDP_vec_target(i);
%     inputs_struc.h = h_hat_vec(i);
% 
%     [Output_struc,inputs_struc] = solve_counterfactual(inputs_struc,fns_struc);
% 
%     % Make sure that errors are around 0
%     d_q_i_C=Output_struc.d_q_i;
%     M_q_C=Output_struc.M_q;
%     pcgdp_prodvec_C_calib_prices = sum(sum(p_q_i_calib.*d_q_i_C.*pdf_A_q_i.*repmat(M_q_C,[1,size(p_q_i_calib,2)])))/N;
%     pcgdp_prodvec_C_rel_calib_prices = pcgdp_prodvec_C_calib_prices/pcgdp_calib_own_prices;
%     error_prod_mult = abs(pcgdp_prodvec_C_rel_calib_prices/pcGDP_vec_target(i) - 1)
%     error_theta = abs(Output_struc.w_S./w_S_target - 1)
% 
%     Output_struc_US_C{i,1} = Output_struc;
%     inputs_struc_US_C{i,1} = inputs_struc;
% 
%     end
% 
%     % Write C across states to excel
%     excel_column = 1;
%     prod_mult = 1;
%     % Write the labels on the right most column
%     xlswrite(file_name,{'poorest','calib','richest'},strcat(sheet_name,';US'),strcat(char('A'+excel_column),'1'));
%     write_label = {'q_1';'h';'prod_mult';'pcGDP';'w_S';'L_U clearing';'L_S_clearing';'ch_sd';'sh_emp_small'};
%     write_label{end+1,1} = 'rel_price';
%     write_label{end+no_q,1} = 'avg_size';
%     write_label{end+no_q,1} = 'theta';
%     write_label{end+no_q,1} = 'rel_prod';
%     write_label{end+no_q,1} = 'sh S entry';
%     xlswrite(file_name,write_label,strcat(sheet_name,';US'),'A2');
%     % Write output
%     xlswrite(file_name,[Output_struc_US_C{1,1}.write_excel,Output_struc_calib.write_excel],strcat(sheet_name,';US'),'B2');
%     diff_model_US = min(Output_struc_US_C{1,1}.size_hist_cdf_model_C) - min(size_hist_cdf_model_calib);
%     sh_exp_less5_US = diff_model_US/(sh_less5_US_data - min(size_hist_cdf_data))
%     xlswrite(file_name,repmat(sh_exp_less5_US,1,size(pcGDP_vec_target,1)+1),strcat(sheet_name,';US'),'B9');
%     xlswrite(file_name,[min(Output_struc_US_C{1,1}.size_hist_cdf_model_C),min(size_hist_cdf_model_calib)],strcat(sheet_name,';US'),'B10');
% end
% 
% 
% end
% end
% 
% 
% % Write calibration prices to stata to check that I am hitting price size
% % slope
% 
% p_q_i_calib_stata = reshape(p_q_i_calib,[no_q*prod_grid_size,1]);
% size_q_i_calib_stata = reshape(size_q_i_calib,[no_q*prod_grid_size,1]);
% pdf_A_q_i_stata = reshape(pdf_A_q_i,[no_q*prod_grid_size,1]);
% M_q_calib_stata = reshape(repmat(M_q_calib,1,prod_grid_size),[no_q*prod_grid_size,1]);
% data_stata = [p_q_i_calib_stata,size_q_i_calib_stata,pdf_A_q_i_stata,M_q_calib_stata];
